#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <CL/cl.h>
#include "fuel_cell_infos.h"
#include "main_matrix.h"
#include "matrix.h"
#include "opencl_tools.h"
#include "params.h"
#include "point.h"
#include "physical_constants.h"
#include "list.h"
#include "export.h"
#include "legend.h"

void FuelCellMain();
void ComputePackedList();
void ComputeStandart();
void ComputeStandartX(int x);
void ComputeDataToExchange();
void Compute(double diff);
void ComputeStep();
void ExchangeResults();
void ExchangeResultsWait();
void RetrieveResults();
void RetrieveResultsMaster();
void RetrieveResultsSlave();
int GetRetrieveResultsTransfertSize(int rank);

void MainMatrixBackup();
void MainMatrixResultSend();
void MainMatrixResultReceive();
void MainMatrixInitWithBackup();
void MainMatrixInitWithBackupBlock(int bx, int bbx, int by, int bby, int bz, int bbz,
        int bsx, int bbsx, int bsy, int bbsy, int bsz, int bbsz);
tData OverlapOneDim(int iPos1, int iPos2, int iSize1, int iSize2);
tData OverlapVolume(Point3 pos1, Point3 size1, Point3 pos2, Point3 size2);
tData NewTemperature(Point3 newPos, Point3 newSize, Point3 oldBlockPos, Point3 oldSize);
void ComputeSymmetric();

int symmetric = 1;
tData ***backupMatrix = NULL;
ListInt *backupBlockNbNodesX = NULL;
ListInt *backupBlockNbNodesY = NULL;
ListInt *backupBlockNbNodesZ = NULL;

tData tMin(tData v1, tData v2){
    if( v1 < v2 )
        return v1;
    return v2;
}

tData tMax(tData v1, tData v2){
    if( v1 > v2 )
        return v1;
    return v2;
}

int main(int argc, char **argv){
/**************************** st Open CL ****************************/
    opencl_init();
    messagecl("Platform Name=%s",clPlatformName);
    double start_time=getTime();
    FuelCellMain();
	messagecl("Total time : %f", getTime() - start_time);

/**************************** en Open CL ****************************/

    return 0;
}

void FuelCellMain(){
    int granularity[3] = {1000, 1000, 1000};

    int i=0;
    do{
    messagecl("Starting computation with granularityX = %d, granularityY = %d, granularityZ = %d",
            granularity[0], granularity[1], granularity[2]);
        InitFuelCellInfos(granularity);
        MainMatrixInit();
        MainMatrixInitWithBackup();
        ExportInit();
        BuildLegendFile();
        double diff = 11. * (i + 1);
        if(granularity[0] == 17)
            diff = 0.1;	//diff = 0.0001; // Temporary hack
// we should change compute function to OpenCL
        Compute(diff);
        MainMatrixBackup();
        granularity[0] += 1;
        granularity[1] += 1;
        granularity[2] += 1;
        i+=3;
    }while(i < 200);

    MainMatrixFree();
}

void MainMatrixBackup(){
    backupMatrix = (tData ***)Matrix3Allocate(zero3, matrixSize, sizeof(tData));

    backupBlockNbNodesX = blockNbNodesX;
    backupBlockNbNodesY = blockNbNodesY;
    backupBlockNbNodesZ = blockNbNodesZ;

    blockNbNodesX = NULL;
    blockNbNodesY = NULL;
    blockNbNodesZ = NULL;
        Matrix3Copy(currentMatrix, backupMatrix, zero3, matrixSize);
        MainMatrixResultSend();
}

void MainMatrixResultSend(){
    int transfertSize = matrixSize.x * matrixSize.y * matrixSize.z;
}

void MainMatrixInitWithBackup(){
    if( backupMatrix == NULL )
        return;
    ListPrint(backupBlockNbNodesX);
    ListPrint(blockNbNodesX);

    int bsx, bbsx;
    int bsy, bbsy;
    int bsz, bbsz;

    int bx = 0;
    int bbx = 0;
    for(int lx=0; lx<blockNbNodesX->size; lx++){
        int by = 0;
        int bby = 0;
        bsx = blockNbNodesX->list[lx];
        bbsx = backupBlockNbNodesX->list[lx];
        for(int ly=0; ly<blockNbNodesY->size; ly++){
            int bz = 0;
            int bbz = 0;
            bsy = blockNbNodesY->list[ly];
            bbsy = backupBlockNbNodesY->list[ly];
            for(int lz=0; lz<blockNbNodesZ->size; lz++){
                bsz = blockNbNodesZ->list[lz];
                bbsz = backupBlockNbNodesZ->list[lz];
                // HERE GO THE CODE
                MainMatrixInitWithBackupBlock(bx, bbx, by, bby, bz, bbz,
                        bsx, bbsx, bsy, bbsy, bsz, bbsz);

                bz += bsz;
                bbz += bbsz;
            }
            by += bsy;
            bby += bbsy;
        }
        bx += bsx;
        bbx += bbsx;
    }
}

void MainMatrixInitWithBackupBlock(int bx, int bbx, int by, int bby, int bz, int bbz,
        int bsx, int bbsx, int bsy, int bbsy, int bsz, int bbsz){
    int x, y, z;

    Point3 newSize = Point3i(bsx, bsy, bsz);
    Point3 oldBlockPos = Point3i(bbx, bby, bbz);
    Point3 oldSize = Point3i(bbsx, bbsy, bbsz);

    for(x = bx; x < bx + bsx; x++){
        if( x < localMatrixOffset.x )
            continue;
        if( x >= localMatrixOffset.x + localMatrixSize.x )
            continue;

        for(y = by; y < by + bsy; y++){
            for(z = bz; z < bz + bsz; z++){
                Point3 newPos = Point3i(x - bx, y - by, z - bz);

                currentMatrix[x][y][z] = NewTemperature(newPos, newSize, oldBlockPos, oldSize);
            }
        }
    }
}

tData NewTemperature(Point3 newPos, Point3 newSize, Point3 oldBlockPos, Point3 oldSize){
    tData temp = 0;
    tData newBlockVolume = 1. /( (tData)newSize.x * (tData)newSize.y * (tData)newSize.z );

    Point3 oldPos;
    for(oldPos.x = 0; oldPos.x<oldSize.x; oldPos.x++){
        for(oldPos.y = 0; oldPos.y<oldSize.y; oldPos.y++){
            for(oldPos.z = 0; oldPos.z<oldSize.z; oldPos.z++){
                tData oldTemp = backupMatrix[oldPos.x + oldBlockPos.x][oldPos.y + oldBlockPos.y][oldPos.z + oldBlockPos.z];
                if( oldTemp == 0)
                    die_cl("oldTemp!");
                tData ov = OverlapVolume(newPos, newSize, oldPos, oldSize);
                temp += oldTemp * ov / newBlockVolume;
            }
        }
    }

    if( temp == 0 )
        die_cl("temp = %f", temp);

    return temp;
}

tData OverlapVolume(Point3 pos1, Point3 size1, Point3 pos2, Point3 size2){
    tData sX = OverlapOneDim(pos1.x, pos2.x, size1.x, size2.x);
    tData sY = OverlapOneDim(pos1.y, pos2.y, size1.y, size2.y);
    tData sZ = OverlapOneDim(pos1.z, pos2.z, size1.z, size2.z);

    return sX * sY * sZ;
}

tData OverlapOneDim(int iPos1, int iPos2, int iSize1, int iSize2){
    tData size1 = 1. / (tData)iSize1;
    tData pos1 = size1 * (tData)iPos1;
    tData pos1R = pos1 + size1;

    tData size2 = 1. / (tData)iSize2;
    tData pos2 = size2 * (tData)iPos2;
    tData pos2R = pos2 + size2;

    if( size1 > 1. )
        die_cl("Size1 > 1");

    if( size2 > 1. )
        die_cl("Size2 > 1");

    tData left = tMax(pos1, pos2);
    tData right = tMin(pos1R, pos2R);

    if( right < left )
        return 0.;

    return right - left;
}

void Compute(double diff){
    int max_it = 500000000;//;
    double start_time = getTime();
    for(int i = 0; i <= max_it; i+=exchangeDiviser){ 
        if( i % 2500 == 0 ){
                ComputeSymmetric();
                double difference = ExportIt(i);
                messagecl("Step: %7d, Temps de calcul : %3.2f s (%3.0f min), diff = %f", i, getTime() - start_time, (getTime() - start_time) / 60, difference);
                if( difference < diff ){
                    messagecl("Finished...");
                    return;
                }
        }
        ComputeStep();
    }
        messagecl("Temps de calcul : %3.2f s", getTime() - start_time);
}

void ComputeSymmetric(){
    if( !symmetric )
        return;

    for(int y=0; y<matrixSize.y / 2; y++)
        for(int x=0; x<matrixSize.x; x++)
            for(int z=0; z<matrixSize.z; z++){
                currentMatrix[x][matrixSize.y - y - 1][z] = currentMatrix[x][y][z];
			}
}

void ComputeStep(){
    for(int i = 1; i < exchangeDiviser; i++){
        Matrix3Swap(&currentMatrix, &oldMatrix);

        ComputeStandart(i);
    }

    Matrix3Swap(&currentMatrix, &oldMatrix);
    ComputeDataToExchange();
    ComputeStandart(2 * exchangeDiviser-1);
}

void ComputeDataToExchange(){
    for(int i = 0; i < exchangeDiviser; i++){
        ComputeStandartX(localMatrixOffset.x + exchangeDiviser + i);
        ComputeStandartX(localMatrixOffset.x + localMatrixSize.x - exchangeDiviser - i - 1);
    }
}

void ComputeStandart(int cover){
    for(int x = localMatrixOffset.x + cover; x < localMatrixOffset.x + localMatrixSize.x - cover; x++)
        ComputeStandartX(x);
}

void ComputeStandartX(int x){
    int px = ePackedListX[x];
    int maxY;

    if(symmetric)
        maxY = (localMatrixOffset.y + localMatrixSize.y + 1) / 2;
    else
        maxY = localMatrixOffset.y + localMatrixSize.y - 1;

    for(int y = localMatrixOffset.y + 1; y < maxY; y++){
        int py = ePackedListY[y];
        for(int z = localMatrixOffset.z + 1; z < localMatrixOffset.z + localMatrixSize.z - 1; z++){
            int pz = ePackedListZ[z];
            CoeffBlock *cb = &blocksCoeff[px][py][pz];
            if( cb->special_type == 0 ){
                /* pas un gaz */
                currentMatrix[x][y][z] =
                    oldMatrix[x][y][z] * cb->self +
                    oldMatrix[x-1][y][z] * cb->transfert[0] +
                    oldMatrix[x+1][y][z] * cb->transfert[1] +
                    oldMatrix[x][y-1][z] * cb->transfert[2] +
                    oldMatrix[x][y+1][z] * cb->transfert[3] +
                    oldMatrix[x][y][z-1] * cb->transfert[4] +
                    oldMatrix[x][y][z+1] * cb->transfert[5] +
                    cb->constant;
                if(currentMatrix[x][y][z] == 0.){
                  int nb_neighbours = 0;
                  if(oldMatrix[x-1][y][z] != tMaterial[MAT_OUTSIDE].Tinit){
                    currentMatrix[x][y][z] += oldMatrix[x-1][y][z];
                    nb_neighbours++;
                  }
                  if(oldMatrix[x+1][y][z] != tMaterial[MAT_OUTSIDE].Tinit){
                    currentMatrix[x][y][z] += oldMatrix[x+1][y][z];
                    nb_neighbours++;
                  }
                  if(oldMatrix[x][y-1][z] != tMaterial[MAT_OUTSIDE].Tinit){
                    currentMatrix[x][y][z] += oldMatrix[x][y-1][z];
                    nb_neighbours++;
                  }
                  if(oldMatrix[x][y+1][z] != tMaterial[MAT_OUTSIDE].Tinit){
                    currentMatrix[x][y][z] += oldMatrix[x][y+1][z];
                    nb_neighbours++;
                  }
                  if(oldMatrix[x][y][z-1] != tMaterial[MAT_OUTSIDE].Tinit){
                    currentMatrix[x][y][z] += oldMatrix[x][y][z-1];
                    nb_neighbours++;
                  }
                  if(oldMatrix[x][y][z+1] != tMaterial[MAT_OUTSIDE].Tinit){
                    currentMatrix[x][y][z] += oldMatrix[x][y][z+1];
                    nb_neighbours++;
                  }
                  if(nb_neighbours)
                    currentMatrix[x][y][z] = currentMatrix[x][y][z] / (tData)nb_neighbours;
                  else
                    currentMatrix[x][y][z] = oldMatrix[x-1][y][z];
                }
                if( ! ( currentMatrix[x][y][z] > 5 && currentMatrix[x][y][z] <= 1000.) ){
                    printf("temp, me: %e,  v1: %e, v2: %e, v3: %e, v4: %e, v5: %e, v6: %e\n",oldMatrix[x][y][z], oldMatrix[x-1][y][z], oldMatrix[x+1][y][z], oldMatrix[x][y-1][z], oldMatrix[x][y+1][z], oldMatrix[x][y][z-1], oldMatrix[x][y][z+1]);
                    die_cl("(not gaz) temperature incorrecte, value: %e, self: %e, Transferts: (%e, %e, %e, %e, %e, %e), constant: %e, coords: (%d, %d, %d)", currentMatrix[x][y][z], cb->self, cb->transfert[0], cb->transfert[1], cb->transfert[2], cb->transfert[3], cb->transfert[4], cb->transfert[5], cb->constant, x, y, z);
                }
            }else{
                /* gaz */
                currentMatrix[x][y][z] = oldMatrix[x][y][z] + (- cb->transfert[0] - cb->transfert[2] - cb->transfert[4]) * oldMatrix[x][y][z] +
                    (
                    oldMatrix[x-1][y][z] * cb->transfert[0] +
                    oldMatrix[x][y-1][z] * cb->transfert[2] +
                    oldMatrix[x][y][z-1] * cb->transfert[4]);
                if( ! ( currentMatrix[x][y][z] > 5 && currentMatrix[x][y][z] < 700) ){
                    printf("temp, me: %e,  v1: %e, v2: %e, v3: %e\n",oldMatrix[x][y][z], oldMatrix[x+1][y][z], oldMatrix[x][y+1][z],  oldMatrix[x][y][z+1]);
                    die_cl("(gaz) temperature incorrecte, value: %e, self: %e, Transferts: (%e, %e, %e), constant: %e, coords: (%d, %d, %d)", currentMatrix[x][y][z], cb->self, cb->transfert[0], cb->transfert[2], cb->transfert[4], cb->constant, x, y, z);
                }
            }
        }
    }

    if( symmetric ){
        maxY--;
        if( localMatrixSize.y % 2 == 1 ){
            int y1 = maxY - 1;
            int y2 = maxY + 1;
            for(int z = localMatrixOffset.z + 1; z < localMatrixOffset.z + localMatrixSize.z - 1; z++)
                currentMatrix[x][y2][z] = currentMatrix[x][y1][z];
        }else{
            int y1 = maxY;
            int y2 = maxY + 1;
            for(int z = localMatrixOffset.z + 1; z < localMatrixOffset.z + localMatrixSize.z - 1; z++)
                currentMatrix[x][y2][z] = currentMatrix[x][y1][z];
        }
    }
}

int GetRetrieveResultsTransfertSize(int rank){
    Point3 lMatrixSize = ComputeLocalMatrixSize(rank);
    int nbExchangeDiviserToIgnore = 1;
    int transfert_size = matrixSize.y * matrixSize.z *
        (lMatrixSize.x - exchangeDiviser * nbExchangeDiviserToIgnore);

    return transfert_size;
}

